% A4_regrid_mm5.m 9/27/2011 Parker MacCready
%
% this takes the processed mm5 files and...
% - interpolates them to the same plaid lat-lon grid
% - stores things in reasonable size chunks (monthly)

% now including wrf files!

clear
tic

% *** USER DEFINED VALUES **********************
dth = 120; % time step for forcing fields (hours) [120, 12, 3, 1]
period_list = {'jun01a';'jun01b';'aug02a';'aug02b'; ...
    'dec03a';'dec03b';'may06a';'may06b';'wrf1a';'wrf1b'};
year_vec = 2001:2009; % span of years to regrid (could be just one)
% year_vec = 2001:2003; % span of years to regrid (could be just one)
% year_vec = 2004:2005; % span of years to regrid (could be just one)
% year_vec = 2006:2007; % span of years to regrid (could be just one)
% year_vec = 2008:2009; % span of years to regrid (could be just one)
% *** END USER DEFINED VALUES ******************

switch dth
    case 120
        indir0 = 'processed_five_day/';
    case 12
        indir0 = 'processed_half_day/';
    case 3
        indir0 = 'processed_three_hour/';
    case 1
        indir0 = 'processed_hourly/';
    otherwise
        disp('that dth is not coded for')
        break
end

% create the output directory
outdir0 = strrep(indir0,'processed','regridded');
if ~(exist(outdir0,'dir')==7); mkdir(outdir0); end

disp(' '); disp(['*** Regridding from ',indir0,' ***']);

% design the plaid grid to interpolate to
lon1 = -131; lon2 = -121.7; lat1 = 40; lat2 = 52;
% aim for 4 km resolution everywhere
RE = 6371e3; % approximate Earth radius (m)
dy = RE*pi*(lat2-lat1)/180;
mlr = pi*(mean([lat1 lat2]))/180;
dx = cos(mlr)*RE*pi*(lon2-lon1)/180;
if 0 % 4 km grid
    ny = round(dy/4000);
    nx = round(dx/4000);
    ilon = linspace(lon1,lon2,nx);
    ilat = linspace(lat1,lat2,ny);
else % stretched grid (has only 60% of the points of the 4 km grid)
    dlon12 = 12e3*180/(RE*cos(mlr)*pi);
    dlon4 = 4e3*180/(RE*cos(mlr)*pi);
    dlat12 = 12e3*180/(RE*pi);
    dlat4 = 4e3*180/(RE*pi);
    iii = 1; % a counter
    ilon = lon1;
    while ilon(iii) <= lon2
        iii = iii + 1;
        ilon = [ilon; ilon(iii-1) + dlon4 + ...
            0.5*(dlon12-dlon4)*(1-tanh(2*(ilon(iii-1)+127.5)))];
    end
    iii = 1; % a counter
    ilat = lat1;
    while ilat(iii) <= lat2
        iii = iii + 1;
        ilat = [ilat; ilat(iii-1) + dlat4 + ...
            0.5*(dlat12-dlat4)* (1+max(-tanh((2*(ilat(iii-1)-41.5))), ...
            (tanh(2*(ilat(iii-1)-50.5)))))];
    end
end
[LON,LAT] = meshgrid(ilon,ilat);
[NR,NC] = size(LON);

% initialize a cell array to hold all the time vectors
np = length(period_list);
td = cell(np,1);
% and arrays for the month and year
yr = cell(np,1);
mo = cell(np,1);
% then get all available time vectors (all periods)
for pp = 1:np
    period = period_list{pp};
    indir = [indir0,period,'/'];
    load([indir,'t2_d2.mat'],'tdvec'); % those from d3 should be identical
    td{pp} = tdvec(:);
    yvec = str2num(datestr(tdvec,'yyyy'));
    yr{pp} = yvec;
    mvec = str2num(datestr(tdvec,'mm'));
    mo{pp} = mvec;
end % end of period loop

% LOOPS

% choose a variable
varlist = {'psfc';'t2';'qair';'rain';'swdown';'lwdown';'u10r';'v10r'};
for vv = 1:length(varlist);
    var = varlist{vv};
    
    % choose a year
    for year = year_vec
        
        outdir = [outdir0,num2str(year),'/'];
        if ~(exist(outdir,'dir')==7); mkdir(outdir); end
        
        % choose a month
        for month = 1:12
            
            disp([' -- working on ',var,' for ',num2str(year),'/', ...
                num2str(month)]);
            
            % find the total number of times that are in
            % the chosen year and month
            nt = [];
            for pp =1:np
                nt(pp) = sum(yr{pp}==year & mo{pp}==month);
            end
            NT = sum(nt);
            
            % create the output matrix
            TD = NaN * ones(NT,1);
            VAR = NaN * ones(NT,NR,NC);
            
            % step through the input files to populate the output matrix
            it = 0; % initialize a counter
            % interpolate the d2 field to the new grid
            for pp = 1:np
                period = period_list{pp};
                ind = find(yr{pp}==year & mo{pp}==month);
                if ~isempty(ind)
                    indir = [indir0,period,'/'];
                    load([indir,var,'_d2.mat']);
                    for ii = 1:length(ind)
                        it = it + 1; % a counter
                        TD(it) = tdvec(ind(ii));
                        oo = squeeze(outmat(ind(ii),:,:));
                        if ~isnan(oo(1,1))
                            F = TriScatteredInterp(lon(:),lat(:),oo(:));
                            ioo = F(LON,LAT);
                            % extrapolate to the north (a tiny triangle)
                            for cc = 1:size(ioo,2)
                                this_col = ioo(:,cc);
                                inan = find(isnan(this_col));
                                if ~isempty(inan);
                                    this_col(inan) = this_col(inan(1)-1);
                                    ioo(:,cc) = this_col;
                                end
                            end
                            
                            VAR(it,:,:) = ioo;
                        end
                    end
                end
            end
            
            % and interpolate in the d3 fields where we have them
            it = 0; % initialize a counter
            for pp = 1:np
                period = period_list{pp};
                ind = find(yr{pp}==year & mo{pp}==month);
                if ~isempty(ind)
                    indir = [indir0,period,'/'];
                    load([indir,var,'_d3.mat']);
                    for ii = 1:length(ind)
                        it = it + 1; % a counter
                        oo = squeeze(outmat(ind(ii),:,:));
                        if ~isnan(oo(1,1))
                            F = TriScatteredInterp(lon(:),lat(:),oo(:));
                            ioo = F(LON,LAT);
                            % fix the weird edge values of the
                            % interpolated field
                            for rr = 1:size(ioo,1)
                                this_row = ioo(rr,:);
                                inan = find(~isnan(this_row));
                                if ~isempty(inan);
                                    ioo(rr,inan(1):(inan(1)+1)) = NaN;
                                end
                            end
                            for cc = 1:size(ioo,2)
                                this_col = ioo(:,cc);
                                inan = find(~isnan(this_col));
                                if ~isempty(inan);
                                    ioo(inan(1):(inan(1)+1),cc) = NaN;
                                end
                            end
                            ivar = squeeze(VAR(it,:,:));
                            mask = ~isnan(ioo);
                            ivar(mask) = ioo(mask);
                            VAR(it,:,:) = ivar;
                            
                        end
                    end
                end
            end
            
            % save the results
            % first create a flag vector to use later when we do the
            % interpolation over missing times
            flag_interpolated = zeros(length(TD),1);
            save([outdir,var,'_',num2str(month),'.mat'], ...
                'VAR','LON','LAT','TD','flag_interpolated','-v7.3');
            
        end % end of month loop
        
    end % end of year loop
    
end % end of variable loop

tictoc = toc;
disp(['  ',num2str(round(tictoc)),' sec to regrid']);





